#!/usr/bin/env python
"""
Vizualize mesh partitionsi using gnuplot or vtk

GNUPLOT:
Build a gnuplot data file and the corresponding gnuplot script (need gnuplot installed)

VTK:
Render the scene on screen (need vtk installed)

Three input files are needed:
.node contains the nodes
.part contains the partitions of the nodes
.edge contains the edges
Write two different files:
- .gpd: contains the data to draw
- .gps: contains the gnuplot script that read the gnuplot data text

USAGE:

vizu.py [options] .node .part .edge
it is very important to give these two or three files. If you dont want edges, just put an empty file here. All the options have to be BEFORE the 3 input files

OPTIONS:

gpl options start with --gpl

-h, --help Print this help
-t Delete temp files (.gpd and .gps). Dont use this option for now, since the automatic plotting does NOT work
-p Do NOT delete temp files
--gpl_out_data arg .gpd file name
--gpl_out_script arg .gps file name
-o, --gpl_output arg pdf output file (drawing) name
-d, --debug Print some debug infos
-i, --in_format arg Input format (for part file in particulat): arg can be triangle or metis (DEFAULT=Metis)
-r, --out_render arg Output renderer: arg can be gnuplot or vtk (DEFAULT=gnuplot)
--vtk_rsphere arg Specifies radius of the nodes (DEFAULT: 0.2)
--vtk_rtube arg Specifies radius of the edges (DEFAUTL: vtk_rsphere/15)
--vtk_filename arg Name of the file to print (DEFUALT: test)
--vtk_pdf, --vtk_png, --vtk_nooutput, --vtk_rasterpdf Format of the output file (DEFAULT: vtk_nooutput). vtk_nooutput means interactive vizualization


FILES FORMAT
.node format
# #_of_nodes dimension (2 or 3)
# node_id x y (z)
9 2 #number of nodes, dimension
1 0 0 0 # some number, x, y, z (if dim=3)
2 1 0 0
3 2 0 0
4 0 1 0
5 1 1 0
6 2 1 0
7 0 2 0
8 1 2 0
9 2 2 0

.part format (Triangle)
# hophophop
9 4 #Number of nodes, number of domains
1 1 # Node, domain
2 1
3 4
4 1
5 1
6 3
7 2
8 2
9 2

.part format (Metis)
#hiphip
1 # domain #. ith line represents node i
2
2
1
4
3
4
4
2

.edge format
#Edges
5 0 #Number of edges, 0
1 1 5 #A number, node1, node2
2 8 9
3 5 8
4 9 6
5 2 3


EXAMPLE:

gnuplot mytest.gps
display mytest.pdf

vtk:
  ./showmesh.py -i metis -r vtk --vtk_rsphere 0.6 --vtk_filename test1 --vtk_pdf ../data/testix.1.node ../data/testix.1.edge ../data/parts_metis/testix.1.mesh.ngraph.part.7

BUGS

Gnuplot: Currently two major problems:
- the 3D drawing works but is ugly. Impossible to get the correct colors and the same kind of point for each node.
- the popen for launching gnuplot does not want to work. Signals an error if the given file does not exist, and thats it, asolutely no more, the .pdf output file is not even created
So, for now, just run \"gnuplot .gps\" after this script finishes to get the image.

Basile Clout, September 2007
"""


import os
import sys
import string
import getopt
import pprint
from subprocess import *
from datetime import *
from string import Template

import vtk

def pnode(fnode, fl_format):
    """ Parse the .node file
    and returns the correspondind
    list of the nodes
    """

    nodes = []
    count = 0
    dim = 0
    nbnodes = 0

    for line in fnode:

        if line[0]=='#':
            continue
        else:
            parts = line.split()

            if count == 0:
                nbnodes = int(parts[0])
                dim = int(parts[1])
                count += 1
                if dim not in [2,3]:
                    print "Error dim in .node file!"

            else:
                if dim == 3:
                    nodes.append(tuple([float(x) for x in parts[1:]]))
                else:
                    nodes.append(tuple([float(x) for x in parts[1:3]]))

    if (len(nodes) != nbnodes):
        print "Error: actual number of nodes does not corresponds " + str(len(nodes)), nbnodes
        sys.exit(2)
                                 
    return nodes, dim


def pedges(fedge, fl_format):
    """ Parse the .edge file and returns the corresponding
    list of edges
    """

    edges = []
    count = 0
    nbedges = 0
    boundary = 0
    for line in fedge:

        if line[0]=='#':
            continue
        else:
            parts = line.split()

            if count == 0:
                nbedges = int(parts[0])
                boundary = int(parts[1])
                count += 1
                if boundary not in (0,1):
                    print "Error boundary in .edge file!"

            else:
                if boundary == 1:
                    edges.append(tuple([int(x) for x in parts[1:]]))
                else:
                    edges.append(tuple([int(x) for x in parts[1:3]]))

    if (len(edges) != nbedges):
        print "Error: actual number of edges does not corresponds " + str(len(edges)), str(nbedges)
        sys.exit(2)

    return edges


def pparts(fparts, fl_format):
    """ Build a list of the partitions:
    [(node, part), (node, part), ...]
    """

    nbnodes = 0
    nbparts = 0
    count = 0
    mapping = []
    domains = []
    i = 0
    p = 0
    dom = 0

    for line in fparts:

        if line[0]=='#':
            continue
        
        else:
            parts = line.split()

            if fl_format == 0: # Check the first line if in triangle format
            
                if count == 0 :
                    nbnodes = int(parts[0])
                    nbparts = int(parts[1])
                    count += 1
                    
                else:
                    if int(parts[1]) not in domains:
                        domains.append(int(parts[1]))
                        
                    mapping.append(tuple([int(x) for x in parts]))
                        
            elif fl_format == 1:

                dom = int(parts[0]) + 1
                if dom not in domains:
                    domains.append(dom)
                    p = p+1

                i = i+1
                mapping.append((i, dom))

    if(fl_format==1):
        nbnodes = i;
        nbparts = p;
                    
    if (len(mapping) != nbnodes ):
        print "Error: actual number of nodes does not corresponds " + str(len(mapping)), nbnodes
        sys.exit(2)
    if (len(domains) != nbparts):
        print "Error: actual number of domains does not corresponds " + str(len(domains)), nbparts
        sys.exit(2)

    domains.sort()
    return mapping, domains


def edgeType(edges, mapping):
    """ Return true if the edge is a cutting edge
    """
    cutedges = []
    oedges = []
    for edge in edges:
        if mapping[edge[0]-1][1] != mapping[edge[1]-1][1]:
            cutedges.append(edge)
        else:
            oedges.append(edge)

    return cutedges, oedges
    

def pcoord(coord, dim):
    """ Returns a string containing the formatting of
    the coordinates from coord, given the dimension (2 or 3)
    """

    if dim == 3:
        return str(coord[0]) + " " + str(coord[1]) + " " + str(coord[2]) +"\n"
    else:
        return str(coord[0]) + " " + str(coord[1]) + "\n"



def build_gpd(nodes, dim, edges, domains, mappings, out):
    """ Build the gnuplot data file
    Returns the string to print in this file
    """

    intro = "# Gnuplot data for mesh vizualization\n"
    intro += "# Generated by " + sys.argv[0] + " on " + str(datetime.now()) + " from the files " + str(sys.argv[-3:]) + "\n"
    intro += "# The Nodes paragraph contains the nodes, the Edges one the ... edges description\n\n"
    
    intro += "# Nodes: id x y z\n"
    out.write(intro)

    for i in domains:
        out.write("# Domain " + str(i) + "\n")
        for (node, domain) in mappings:
            if i==domain:
                out.write(pcoord(nodes[node-1], dim))

        out.write("\n\n")
        
    
    out.write("# Edges: node1 node2\n")

    # Cutedges list contains edges between two nodes within two different domains
    cutedges, oedges = edgeType(edges, mappings) 
    #print "cutedges:"
    #pprint.pprint(cutedges)
    #print "other edges:"
    #pprint.pprint(oedges)
    
    intro += "# No cutting-edges\n"
    for edge in oedges:
        idnode1 = nodes[edge[0]-1]
        idnode2 = nodes[edge[1]-1]
        out.write(pcoord(idnode1, dim))
        out.write(pcoord(idnode2, dim))
        out.write("\n")

    out.write("\n")
    out.write("# Cutting edges\n")

    for edge in cutedges:
        idnode1 = nodes[edge[0]-1]
        idnode2 = nodes[edge[1]-1]
        out.write(pcoord(idnode1, dim))
        out.write(pcoord(idnode2, dim))
        out.write("\n" )


def build_gps(nodes, dim, edges, domains, mappings, filename, out_gpd, out):
    """ Build the gnuplot script file
    Returns the string to print in this file
    """
    
    intro = "# Gnuplot script for mesh vizualization\n"
    intro += "# Generated by " + sys.argv[0] + " on " + str(datetime.now()) + " from the files " + str(sys.argv[-3:]) + "\n"
    intro += "set term pdf\n"
    intro += "set output \"%s\"\n" % filename
    intro += "set key right bottom\n"
    intro += "set title \'Partition of the mesh in " + str(len(domains)) + " subdomains\'\n"
    intro += "set pointsize 2.0\n"
    xmin, xmax, ymin, ymax = window(nodes)
    intro += "set xrange [%d:%d]\n" % (xmin, xmax)
    intro += "set yrange [%d:%d]\n" % (ymin, ymax)
    
    out.write(intro)

    out.write("\n# Plot\n")
    if dim==3:
        plot = "splot "
        myrange = "1:2:3"
    else:
        plot = "plot "
        myrange = "1:2"
        
    out.write(plot)

    tnodes = Template("\'" + out_gpd +"\' index $index using $urange title '$title' with points pt 1 lw 3 lc $color$end")
    tedges= Template("\'" + out_gpd +"\' index $index using $urange title '$title' with l lw $lwsize lc $color$end")
    
    # Print nodes
    for i in domains:
        out.write(tnodes.substitute(index = i-1, color = i, urange = myrange, title = "domain " + str(i), end = ",\\\n"))

    # Print edges
    out.write(tedges.substitute(index = i, color = i+1, urange = myrange, lwsize = 1, title = "domain edges", end = ",\\\n"))
    out.write(tedges.substitute(index = i+1, color = i+2, urange = myrange, lwsize = 5, title = "cutting edges", end = ""))

    out.write("\n\n print \"" + filename + " created!\"\n");



def build_vtk(_nodes, _dim, _edges, _domains, _mappings, _rsphere, _rtube, filename, fl_file, fl_legend):
    """ Build vtk rendering from standard objects
    """
    
    xmin, xmax, ymin, ymax = window(_nodes)
    wx, wy = 800, 500                   # Rendering window size (and size of the printed graphics)
    colors = [   
               [0.2, 0.2, 1.0, 1.0],    # Deep Blue
               [1.0, 0.0, 0.0, 1.0],    # Red
               [1.0, 1.0, 0.0, 1.0],    # Yellow
               [0.0, 1.0, 0.0, 1.0],    # Green
               [1.0, 0.5, 0.0, 1.0],    # Orange
               [1.0, 0.0, 1.0, 1.0],    # Magenta
               [0.5, 0.5, 0.5, 1.0],    # Grey
               [0.0, 0.0, 0.0, 1.0],    # Black
               [0.0, 1.0, 1.0, 1.0],    # Cyan
               [0.0, 0.5, 0.5, 1.0],    # Turquoise
               [1.0, 1.0, 1.0, 1.0],    # White
               [0.5, 1.0, 0.5, 1.0],    # Light green
             ]


  
    lut = vtk.vtkLookupTable()
    lut.SetNumberOfColors(len(colors))
    lut.Build()
    idx=0
    for i in colors:
        lut.SetTableValue(idx, i)
        idx = idx+1 
    
    # Nodes (Points)
    nodes = vtk.vtkPoints()
    nodes.SetDataTypeToFloat()
    nodesColor = vtk.vtkIdTypeArray()
    nodesCell = vtk.vtkCellArray()
    i=0
    
    for node in _nodes:

        if(_dim==2):
            id = nodes.InsertPoint(i, node[0], node[1], 0.0)
        else:
            id = nodes.InsertPoint(i, node[0], node[1], node[2])
        
        nodesColor.InsertNextValue(_mappings[i][1])    # Color (domain) corresponding to this node

        nodesCell.InsertNextCell(1)
        nodesCell.InsertCellPoint(i)
        i = i+1

    # Add nodes corresponding to the legend spheres
    if(fl_legend!=0):
        
        for part in _domains:
            id = nodes.InsertPoint(i, xmax+(xmax-xmin)/10, ymin+(part-1+2)*(ymax-ymin)/10, 0.0)
            nodesColor.InsertNextValue(part)
            nodesCell.InsertNextCell(1)
            nodesCell.InsertCellPoint(i)
            i = i+1
        
    nodesPoly = vtk.vtkPolyData()
    nodesPoly.SetPoints(nodes)
    nodesPoly.GetPointData().SetScalars(nodesColor)
   # nodesPoly.SetVerts(nodesCell)

    sphere = vtk.vtkSphereSource()
    sphere.SetThetaResolution(16)
    sphere.SetPhiResolution(32)
    sphere.SetRadius(_rsphere)

    spheres = vtk.vtkGlyph3D()
    spheres.SetInput(nodesPoly)
    for i in _domains:
        spheres.SetSource(i-1, sphere.GetOutput())
    spheres.SetColorModeToColorByScalar()
    spheres.SetScaleModeToDataScalingOff()
#    spheres.SetIndexModeToScalar()
    spheres.SetRange(0, len(colors))

    #Edges
    cutedgeColors = vtk.vtkIdTypeArray()
    cutedgeArray = vtk.vtkCellArray()
    cutedgePoints = vtk.vtkPoints()
    cutedgePoints.SetDataTypeToFloat()

    cutedges, oedges = edgeType(_edges, _mappings)
    
    for edge in cutedges:
        if _dim == 2 :
            pt1 = cutedgePoints.InsertNextPoint(_nodes[edge[0]-1][0], _nodes[edge[0]-1][1], 0.0)
            pt2 = cutedgePoints.InsertNextPoint(_nodes[edge[1]-1][0], _nodes[edge[1]-1][1], 0.0)
        else:
            pt1 = cutedgePoints.InsertNextPoint(_nodes[edge[0]-1])
            pt2 = cutedgePoints.InsertNextPoint(_nodes[edge[1]-1])

        cutedgeColors.InsertNextValue(2)

        cutedgeArray.InsertNextCell(2)
        cutedgeArray.InsertCellPoint(pt1)
        cutedgeArray.InsertCellPoint(pt2)


    if fl_legend != 0:
        pt1 = cutedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10, ymin, 0.0)
        pt2 = cutedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10+(xmax-xmin)/22, ymin, 0.0) 
        cutedgeColors.InsertNextValue(2)
        cutedgeArray.InsertNextCell(2)
        cutedgeArray.InsertCellPoint(pt1)
        cutedgeArray.InsertCellPoint(pt2)
        
    cutedgesPolyData = vtk.vtkPolyData()
    cutedgesPolyData.SetPoints(cutedgePoints)
    cutedgesPolyData.SetLines(cutedgeArray)
    cutedgesPolyData.GetCellData().SetScalars(cutedgeColors)


    # Other edges
    oedgeColors = vtk.vtkIdTypeArray()
    oedgeArray = vtk.vtkCellArray()
    oedgePoints = vtk.vtkPoints()
    oedgePoints.SetDataTypeToFloat()

    if fl_file != 3:
        oecolor = 10
    else:
        oecolor = 1
    
    for edge in oedges:
        if _dim == 2 :
            pt1 = oedgePoints.InsertNextPoint(_nodes[edge[0]-1][0], _nodes[edge[0]-1][1], 0.0)
            pt2 = oedgePoints.InsertNextPoint(_nodes[edge[1]-1][0], _nodes[edge[1]-1][1], 0.0)
        else:
            pt1 = oedgePoints.InsertNextPoint(_nodes[edge[0]-1])
            pt2 = oedgePoints.InsertNextPoint(_nodes[edge[1]-1])

        oedgeColors.InsertNextValue(oecolor)
        
        oedgeArray.InsertNextCell(2)
        oedgeArray.InsertCellPoint(pt1)
        oedgeArray.InsertCellPoint(pt2)

    if fl_legend != 0:
        pt1 = oedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10, ymin+(ymax-ymin)/10, 0.0)
        pt2 = oedgePoints.InsertNextPoint(xmax+(xmax-xmin)/10+(xmax-xmin)/22, ymin+(ymax-ymin)/10, 0.0) 
        oedgeColors.InsertNextValue(oecolor)
        oedgeArray.InsertNextCell(2)
        oedgeArray.InsertCellPoint(pt1)
        oedgeArray.InsertCellPoint(pt2)
        
    oedgesPolyData = vtk.vtkPolyData()
    oedgesPolyData.SetPoints(oedgePoints)
    oedgesPolyData.SetLines(oedgeArray)
    oedgesPolyData.GetCellData().SetScalars(oedgeColors)

    tubes = vtk.vtkTubeFilter()
    tubes.SetInput(cutedgesPolyData)
    tubes.SetRadius(5*_rtube)
    tubes.CappingOff()
    tubes.SetNumberOfSides(12)
    tubes.SetGenerateTCoordsToOff()

    otubes = vtk.vtkTubeFilter()
    otubes.SetInput(oedgesPolyData)
    otubes.SetRadius(_rtube)
    otubes.CappingOff()
    otubes.SetNumberOfSides(12)
    otubes.SetGenerateTCoordsToOff()
        
    spheresMap = vtk.vtkPolyDataMapper()
    spheresMap.SetInputConnection(spheres.GetOutputPort())
    spheresMap.SetLookupTable(lut)
    spheresMap.SetScalarRange(1, len(_domains))

    cutedgesMap = vtk.vtkPolyDataMapper()
    cutedgesMap.SetInputConnection(tubes.GetOutputPort())
    cutedgesMap.SetLookupTable(lut)
    cutedgesMap.SetScalarRange(1, 11)
    cutedgesMap.ScalarVisibilityOn()
    oedgesMap = vtk.vtkPolyDataMapper()
    oedgesMap.SetInputConnection(otubes.GetOutputPort())
    oedgesMap.SetLookupTable(lut)
    oedgesMap.SetScalarRange(1, 11)
    oedgesMap.ScalarVisibilityOn()
    
    spheresActor = vtk.vtkActor()
    spheresActor.SetMapper(spheresMap)
    spheresActor.GetProperty().SetSpecularColor(0,1,0)

    cutedgesActor = vtk.vtkActor()
    cutedgesActor.SetMapper(cutedgesMap)
    cutedgesActor.GetProperty().SetOpacity(1.0)
    cutedgesActor.GetProperty().SetInterpolationToGouraud()
    oedgesActor = vtk.vtkActor()
    oedgesActor.SetMapper(oedgesMap)
    oedgesActor.GetProperty().SetOpacity(1.0)
    oedgesActor.GetProperty().SetInterpolationToGouraud()

    #  camera = vtk.vtkCamera()
    #  camera.SetFocalPoint((xmax-xmin)/2 + (xmax-xmin)/3, (ymax-ymin)/2, 0)
    #  camera.SetPosition((xmax-xmin)/2 + (xmax-xmin)/2, (ymax-ymin)/2, 20)
        

    ren = vtk.vtkRenderer()
   # ren.AddActor(nodesActor)
    ren.AddActor(spheresActor)
    ren.AddActor(cutedgesActor)
    ren.AddActor(oedgesActor)

    # Add legend for the domain colors
    for part in _domains:
        text = vtk.vtkVectorText()
        line = "domain "+str(part)
        text.SetText(line)
        mapText=vtk.vtkPolyDataMapper()
        mapText.SetInputConnection(text.GetOutputPort())
        textActor=vtk.vtkFollower()
        textActor.SetMapper(mapText)
        textActor.SetScale(1.5,1.5,1.5)
        textActor.AddPosition(xmax+(xmax-xmin)/10+(xmax-xmin)/20,  ymin+(part-1+2)*(ymax-ymin)/10, 0.0)
        if (fl_file == 3):
            textActor.GetProperty().SetColor(0.0, 0.0, 0.0)
        ren.AddActor2D(textActor)
        ren.ResetCameraClippingRange()
        textActor.SetCamera(ren.GetActiveCamera())

    # Legend for the cutting edges 
    text = vtk.vtkVectorText()
    line = "cutting e. "#(" + str(len(cutedges)) +")"
    text.SetText(line)
    mapText=vtk.vtkPolyDataMapper()
    mapText.SetInputConnection(text.GetOutputPort())
    textActor=vtk.vtkFollower()
    textActor.SetMapper(mapText)
    textActor.SetScale(1.5,1.5,1.5)
    textActor.AddPosition(xmax+(xmax-xmin)/10+(xmax-xmin)/20, ymin, 0.0)
    if (fl_file == 3):
        textActor.GetProperty().SetColor(0.0, 0.0, 0.0)
    ren.AddActor2D(textActor)
    ren.ResetCameraClippingRange()
    textActor.SetCamera(ren.GetActiveCamera())

    # Legend for the other edges
    text = vtk.vtkVectorText()
    line = "other e. "#(" + str(len(oedges)) + ")"
    text.SetText(line)
    mapText=vtk.vtkPolyDataMapper()
    mapText.SetInputConnection(text.GetOutputPort())
    textActor=vtk.vtkFollower()
    textActor.SetMapper(mapText)
    textActor.SetScale(1.5,1.5,1.5)
    textActor.AddPosition(xmax+(xmax-xmin)/10+(xmax-xmin)/20, ymin+(ymax-ymin)/10, 0.0)
    if (fl_file == 3):
        textActor.GetProperty().SetColor(0.0, 0.0, 0.0)
    ren.AddActor2D(textActor)
    ren.ResetCameraClippingRange()
    textActor.SetCamera(ren.GetActiveCamera())


    if fl_file != 3:
	ren.SetBackground(0.1, 0.2, 0.4)
    else:
	ren.SetBackground(1.0, 1.0, 1.0)

    ren.ResetCamera(xmin, xmax, ymin, ymax, 0.0, 0.0)
    
    ren.ResetCameraClippingRange()
    
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(ren)
#    renWin.SetSize(wx, wy)
    renWin.FullScreenOn()

    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    style = vtk.vtkInteractorStyleTrackballCamera()
    iren.SetInteractorStyle(style)
    iren.Initialize()

    


    if (fl_file == 0):                  # Interactive 3D vizualization
        iren.Start()
    elif (fl_file == 1):                # pdf file
        print "Starting normal pdf file generation ..."
        exp = vtk.vtkGL2PSExporter()
        exp.SetRenderWindow(renWin)
        exp.SetFileFormatToPDF()
        #exp.CompressOff()
        #exp.PS3ShadingOff()
        exp.SetFilePrefix(filename)
        exp.Write()
        print filename + ".pdf written!"

    elif (fl_file==2):
        print "Starting raster pdf file generation ..."
        renWin.Render()
        exp = vtk.vtkGL2PSExporter()
        exp.SetRenderWindow(renWin)
        exp.SetFileFormatToPDF()
        exp.SetFilePrefix(filename)
        exp.Write3DPropsAsRasterImageOn()
        exp.Write()
        print filename + ".pdf written!"
                          
    elif (fl_file == 3):                # png file
        print "Starting png file generation ..."
        renWin.Render()
        w2i = vtk.vtkWindowToImageFilter()
        w2i.SetInput(renWin)
        
        writer = vtk.vtkPNGWriter()
        writer.SetInput(w2i.GetOutput())
        writer.SetFileName(filename + ".png")
        writer.Write()
        print filename + ".png written!"
    else:
        print "ERROR: Bad file flag. Exit!"
        sys.exit(1)
        

def window(nodes):
    """ Return the window coordinates:
    in order, xmin, xmax, ymin, ymax
    """

    xmin = 100000000
    xmax = -100000000
    ymin = 100000000
    ymax = -100000000
    for node in nodes:
        if node[0]<xmin:
            xmin = node[0]
        elif node[0]>xmax:
            xmax = node[0]
        if node[1]<ymin:
            ymin = node[1]
        elif node[1]>ymax:
            ymax = node[1]

    return xmin-1, xmax+1, ymin-1, ymax+1
        
        

def main():
    """ Main function
    """
    gnuplot = "/opt/local/bin/gnuplot"
    temp = 0                            # Delete gnuplot data and script files?
    debug = 0                           # Debug flag
    fl_format = 1       # Format of the part file: 0=triangle, 1=metis
    fl_render = 0                       # Render option: 0=gnuplot, 1=vtk

    rsphere = 0.2
    rtube = 0.01
    fl_rtube = 0
    vtk_fl_legend = 1
    
    vtk_filename = "test"                   # vtk filename
    vtk_fl_file = 0                         # vtk type of file: 0=No file (3D interactive vizualization), 1=normal pdf file, 2=raster pdf file, 3=png file
    
    in_node = "test.node"
    in_edge = "test.edge"
    in_part = "test.part"

    out_gpd = "test.gpd"
    out_gps = "test.gps"

    mesh = "test.pdf"
    
    try:
        opts, arg = getopt.getopt(sys.argv[1:], "pfo:dhi:r:", ["help","debug", "in_format=", "out_render", "out_data=", "out_script=", "output=", "vtk_rsphere=", "vtk_rtube=", "vtk_filename=", "vtk_nooutput", "vtk_png", "vtk_pdf", "vtk_rasterpdf", "vtk_legend="])
    except getopt.error, msg:
        print msg
        print "For help, use --help"
        sys.exit(2)

    # Process options
    for o, a in opts:

        if o in ("-h", "--help"):
            print __doc__
            sys.exit(0)
        elif o == "-t":
            temp = 1
        elif o=="-p":
            temp=0
        elif o == "--out_data":
            out_gpd = a
            fl_render = 0
        elif o == "--out_script":
            out_gps = a
            fl_render = 0
        elif o in ("-o", "--output"):
            mesh = a
            fl_render = 0
        elif o in ("-d", "--debug"):
            debug = 1
        elif o == "--vtk_rsphere":
            rsphere = float(a)
            fl_render = 1
        elif o == "--vtk_rtube":
            rtube = float(a)
            fl_rtube = 1
            fl_render = 1
        elif o == "--vtk_filename":
            vtk_filename = a
            fl_render = 1
        elif o == "--vtk_nooutput":
            vtk_fl_file = 0
            fl_render = 1
        elif o == "--vtk_png":
            vtk_fl_file = 3
            fl_render = 1
        elif o == "--vtk_pdf":
            vtk_fl_file = 1
            fl_render = 1
        elif o == "--vtk_rasterpdf":
            vtk_fl_file = 2
            fl_render = 1
        elif o == "--vtk_legend":
            vtk_fl_legend = int(a)
        elif o in ("-i", "--in_format"):
            if a=="triangle":
                fl_format = 0
            elif a=="metis":
                fl_format = 1
            else:
                print "Input format \""+a+"\" unknown!"
                sys.exit(1);
        elif o in ("-r","--out_render"):
            if a=="gnuplot":
                fl_render = 0
            elif a=="vtk":
                fl_render = 1
            else:
                print "Ouput renderer \""+a+"\" unknown!"
                sys.exit(1);
        else:
            print "Option " + o + " unknown. bye ;)"
            sys.exit(2)

    if (fl_rtube != 1):
        rtube = rsphere/15


    if len(arg)<3:
        print "Problem with the given arguments. See help below"
        print __doc__
        sys.exit(2)
    else:
        in_node = arg[0]
        in_edge = arg[1]
        in_part = arg[2]

    try:
        fnode = file(in_node, "r")
        fedge = file(in_edge, "r")
        fpart = file(in_part, "r")
    except IOError:
        print "Impossible to open one of these files: ", in_node, in_edge, in_part
        sys.exit(2)

    print "Use input files :", in_node, in_edge, in_part

    # Build standard objects
    nodes, dim = pnode(fnode, fl_format)
    edges = pedges(fedge, fl_format)
    mappings, domains = pparts(fpart, fl_format)

    if debug == 1:
        print "nodes:"
        pprint.pprint(nodes)
        print "edges"
        pprint.pprint(edges)
        print "mappings"
        pprint.pprint(mappings)
        print "domains"
        pprint.pprint(domains)

    # Gnuplot rendering
    if fl_render == 0:
        ogpd = file(out_gpd, "w")
        ogps = file(out_gps, "w")
        build_gpd(nodes, dim, edges, domains, mappings, ogpd)
        build_gps(nodes, dim, edges, domains, mappings, mesh, out_gpd, ogps)

        order = gnuplot + " " + out_gps
        #p = Popen(order, shell=True, stdout=PIPE, stderr=PIPE)
        #p.wait()
        #print p.stdout.read(), p.stderr.read()
    
        #stat = os.system("/opt/local/bin/gnuplot -persist test.gps")
        print "Operation seems OK! \nScript written in " + out_gps +". Run it with \"gnuplot " + out_gps + "\"."
    
        if temp == 1:
            os.remove(out_gpd)
            os.remove(out_gps)

    # VTK rendering
    elif fl_render == 1:
        build_vtk(nodes, dim, edges, domains, mappings, rsphere, rtube, vtk_filename, vtk_fl_file, vtk_fl_legend)

if __name__ == '__main__':
    main()
    

                                 
    

